AI identifies potent inducers of breast cancer stem cell differentiation based on adversarial learning from gene expression data

Abstract Cancer stem cells (CSCs) are a subpopulation of cancer cells within tumors that exhibit stem-like properties and represent a potentially effective therapeutic target toward long-term remission by means of differentiation induction. By leveraging an artificial intelligence approach solely based on transcriptomics data, this study scored a large library of small molecules based on their predicted ability to induce differentiation in stem-like cells. In particular, a deep neural network model was trained using publicly available single-cell RNA-Seq data obtained from untreated human-induced pluripotent stem cells at various differentiation stages and subsequently utilized to screen drug-induced gene expression profiles from the Library of Integrated Network-based Cellular Signatures (LINCS) database. The challenge of adapting such different data domains was tackled by devising an adversarial learning approach that was able to effectively identify and remove domain-specific bias during the training phase. Experimental validation in MDA-MB-231 and MCF7 cells demonstrated the efficacy of five out of six tested molecules among those scored highest by the model. In particular, the efficacy of triptolide, OTS-167, quinacrine, granisetron and A-443654 offer a potential avenue for targeted therapies against breast CSCs.

Putting the three modules together, the DREDDA architecture aims to solve the so-called unsupervised domain adaptation task (1).Indeed, while ground truth labels (i.e., the cell cluster labels) are only available for the source domain to train the model, the actual aim is to apply the trained model to the target domain (i.e., the LINCS dataset).Compared with Ganin et al. (1), the DREDDA model includes several modifications to make it suitable for gene expression data.
In further detail, to train DREDDA on the unsupervised domain adaptation task, three objective (loss) functions were used.The first objective is the task classification loss for the classification task, which is a cross-entropy between prediction ( ̂) ground truth label (, one-hot encoded): where  ̂() is the output vector of the task classifier   ( () ′ ;   ) and  ̂() takes its -th dimension.The second objective is the adversarial domain loss for the adversarial domain classifier,   () = − () log  ̂() −(1 −  () ) log(1 −  ̂() ) (2) where  () = 0 if the -th example belongs to the source domain and  () The sign in front of   is negative to maximize the domain adversarial classifier loss and thus discourage the use of domain-specific information during the training phase.

Model Training and Implementation
Training of the network using the objective function described in Equation ( 5) can be done with a two-part update per training step.Specifically, one step updates the parameters that minimize the objective ( _ ,  _ ,   ,   ), while the other one updates the parameters that maximize the objective (  ).
Note that gradient ascend is performed on   while gradient descent is performed on all other parameters.The model is implemented using PyTorch 1.8 (4) deep learning framework and can be run on any NVIDIA CUDA-capable GPUs with ≥ 10 memory.

Model Testing and Predictions
After each training epoch, the model's performance was evaluated in the source domain based on multi-class classification loss and accuracy.At the same time, the adversarial classifier was evaluated on the full dataset based on the ability to discriminate between data points coming from the source domain and the target domain.The model was optimized towards the highest source domain accuracy, with target domain accuracy within the 40%-60% range.The trained model was then applied to score each profile in the LINCS database and a final unique drug-related score was obtained by averaging over all the scores obtained for the same drug across different cell lines.

Gene Set Enrichment Analysis of DREDDA-prioritized drugs
Drug-Set Enrichment Analysis (DSEA) (5), which is built on the classical Gene Set Enrichment Analysis (GSEA) (6), identifies pathways consistently dysregulated across a set of experimental conditions (such as drug treatments).DSEA was performed using the top 30 drugs prioritized by DREDDA as the foreground set and the bottom 30 as the background set on the LINCS level 5 dataset for the "C5 GO Biological Process" and "C5 GO Molecular Function" categories as obtained from the MsigDB v7.2 (7).The analysis was performed using the gep2pep R package (8).Summaries for the positive and negative expression within families of pathways were obtained by counting the pathways with positive or negative enrichment scores that were found significant (p < 0.05) by the DSEA analysis and fell below a given level of the GO category.In this case, negative enrichments for GO terms starting with "Negative" (such as "Negative regulation of cell differentiation") were counted as positive.
The drugs, IC50, and working concentrations for cell treatments are shown in Table S5.All drugs were dissolved in DMSO (stock solution) and further diluted in a cell medium (working solution).The negative control was the solvent alone at its final concentration.

Cell viability
Cell viability was assessed by CellTiter 96 AQueous One Solution Cell Proliferation Assay (Promega BioSciences Inc., San Luis Obispo, CA), according to the manufacturer's instructions, 72 hours after cell treatment (5 × 10 3 cells/well, 96-well plates).Each drug concentration was tested in at least three independent experiments for each molecule and each cell line.

Mammosphere assay and CSC self-renewal
Cells were seeded at sub-confluence in 6-well plates.After 16 hrs they were treated with two different concentrations of each compound or DMSO (drug solvent) as negative control.At 24 hrs post-treatment, cells were washed out of the serum and seeded (3000 cells/well) in 24-well low attachment plates with stem cell medium, consisting of 1% B27 (Invitrogen, Carlsbad, CA), 10ng/ml bFGF (Invitrogen), 20 ng/ml EGF (Sigma-Aldrich) in DMEM/F12 (SIGMA-Aldrich) supplemented with 1% penicillin/streptomycin) to form primary (P0) mammospheres.The number and diameter of P0 mammospheres were counted after a period of 7-14 days to estimate the number of treatmentresistant CSCs.Moreover, their diameter was measured to analyze their capacity to proliferate.Then, spheres were disaggregated, and single cells were replated (3000 cells/well) to form secondary (P1) mammospheres.P1 mammospheres were analyzed after further 7-14 days to calculate the CSC self-renewal (P1 mammospheres/P0 plated single cells x 100).

FACS analysis
Adherent cells were treated for 24 hrs with two different concentrations of each compound or negative control (DMSO), then cells were washed with phosphate-buffered saline (PBS) and harvested with 0.05% trypsin/0.025%EDTA.Detached cells were washed and resuspended in PBS supplemented with 2% FBS, 0.2% sodium azide (Staining buffer).Combinations of fluorochromeconjugated monoclonal antibodies against human CD44 (PE; Santa Cruz cat.# SC-18849-PE) and CD24 (FITC; Beckton Dickinson cat.# BD555427) were added to the cell suspension (1 x 10 6 cells) and incubated at 4°C in the dark for 30 min.The labeled cells were washed in a staining buffer and then analyzed on a BD Accuri Flow Cytometer (BD).

Statistics
Ordinary one-way analysis of variance (ANOVA) followed by Tukey's multiple comparisons test, through GraphPad Prism 6 software, was applied for the comparison of groups of experimental data.
Values analyzed are the average +/-SD of at least three independent experiments.

Figure S1
Comparison between DREDDA scores (predicting differentiation) and DECCODE scores (predicting stemness).Although at low DREDDA scores, there is no clear correlation between the two measures, high DREDDA scores tend to correspond to low DECCODE scores and high DECCODE scores tend to correspond to low DREDDA scores.

Figure S2
The number and proportion of genes that are commonly up-(down-) regulated by the top-30 drugs.

Figure S3
Interaction network between the tested drugs and their known targets (analysis performed using the STITCH tool ( 9))     Regulation of systemic arterial blood pressure 0.9 5.79E-08 Potassium ion transport 0.9 5.79E-08 Excretion 0.9 5.79E-08 Regulation of blood pressure 0.9 5.79E-08 Specification of symmetry 0.9 5.79E-08 Regulation of hormone levels 0.9 5.79E-08 Response to auditory stimulus 0.9 5.79E-08 Organic anion transport 0.9 5.79E-08 Organic hydroxy compound transport 0.9 5.79E-08 = 1 if it belongs to the target domain. ̂() is the output of the adversarial domain classifier   ( () ′ ;   ).Inspired by the Deep Domain Confusion (DDC) framework (2), we further introduced a third objective, the domain confusion loss, which explicitly enforces the similarity of intermediate network values between source domain and target domain examples by minimizing a Maximum Mean Discrepancy (MMD) (3) term, where (  ) is the intermediate value in the network after it takes as input   .We chose (  ) to represent the -th intermediate layer of the task classifier, i.e., (  ) =  , (  ′;   ).Only the examples in one mini-batch are used to estimate MMD.The square of this term is used in the final objective function,  =   −   +   = {}, {}) 2

Figure S4
Figure S4 Percentage of cell viability and IC50 of MCF7 and MDA-MB-231 cells after exposure to increasing drug concentrations.MTT test results are reported as a percentage of cell viability.Each value was normalized with respect to its control.Dose-response curves were used to generate IC50.Statistics and IC50 curves were performed using Prism GraphPad 6.0.*, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001 The dependency of each loss term w.r.t. the model parameters can be inferred from the DREDDA architecture.Specifically, they are   =   ( _ ,   ,   ),   =   ( _ ,  _ ,   ,   ) and   =   ( _ ,  _ ,   ,   ).The parameters are optimized via the following minimax objective,   ( _ ,   ,   ) −   ( _ ,  _ ,   ,   ) +   ( _ ,  _ ,   ,   ) For each training step, an equal number of source domain and target domain examples are sampled from the source domain and target domain datasets.Using the sampled examples,   is computed using the sampled source domain examples, while   and   are computed using the sampled source and target domain examples.If plain minibatch stochastic gradient optimization is used, the parameters are updated as follows,  _ ←  _ −   _  _ ←  _ −   _    ←   −       ←   −       ←   +

Table S2
Genes that are commonly up-or down-regulated by the top 30 molecules.

Table S3
Pathways enriched by genes commonly dysregulated by the top 10 drugs.

Table S4
Top-enriched GO terms under the Biological Process category resulting from drug set enrichment analysis (DSEA) of the top 30 drugs vs bottom 30 drugs.

Table S5
Top-enriched GO terms under the Cellular Component category resulting from drug set enrichment analysis (DSEA) of the top 30 drugs vs bottom 30 drugs.